/***************************************************************************
 *   File DFT.cpp		author belot nicolas (karrakis)                    *
 *   define DFT
 *   Manipulation of DFT		                                   *
 ***************************************************************************
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 *   This program is distributed in the hope that it will be useful,       *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
 *   GNU General Public License for more details.                          *
 *                                                                         *
 *   You should have received a copy of the GNU General Public License     *
 *   along with this program; if not, write to the                         *
 *   Free Software Foundation, Inc.,                                       *
 *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
 ***************************************************************************/

#include "DFT.h"
#define max(a,b) a>b?a:b
#define min(a,b) a<b?a:b
namespace libtrckr {
DFT::DFT():size(1024){
	dft=(fftw_complex*)fftw_malloc(sizeof(fftw_complex)*size);
	for(unsigned int i =0;i<size;i++){
		dft[i][0]=0.0;
		dft[i][1]=0.0;
	}
}

DFT::DFT(unsigned int Size):size(Size){
	dft=(fftw_complex*)fftw_malloc(sizeof(fftw_complex)*size);
	for(unsigned int i =0;i<size;i++){
		dft[i][0]=0.0;
		dft[i][1]=0.0;
	}
}

DFT::DFT(double* real,double* img, unsigned int Size):size(Size){
	dft=(fftw_complex*)fftw_malloc(sizeof(fftw_complex)*size);
	for(unsigned int i=0;i<size;i++){
		dft[i][0]=real[i];
		dft[i][1]=img[i];
	}
}

DFT::DFT(fftw_complex* data,unsigned int Size):size(Size){
	dft=(fftw_complex*)fftw_malloc(sizeof(fftw_complex)*size);
	memcpy(dft, data, sizeof(fftw_complex)*size);
	/*for(unsigned int i=0;i<size;i++){
		dft[i][0]=data[i][0];
		dft[i][1]=data[i][1];
	}*/
}

DFT::DFT(const DFT &d2):size(d2.size){
	dft=(fftw_complex*)fftw_malloc(sizeof(fftw_complex)*size);
	memcpy(dft,d2.dft, sizeof(fftw_complex)*size);
	/*for(unsigned int i=0;i<size;i++){
		dft[i][0]=d2.dft[i][0];
		dft[i][1]=d2.dft[i][1];
	}*/
}

DFT::~DFT(){
	fftw_free(dft);
}

unsigned int DFT::getSize() const{
	return size;
}

double *DFT::getRealParts() const{
	double *ret=(double*)malloc(sizeof(double)*size);
	for(unsigned int i=0; i<size;i++){
		ret[i]=dft[i][0];
	}
	return ret;
}

uint DFT::getRealParts(double* data, uint dsize) const{
	uint minsize=min(size,dsize);
	for(uint i=0;i<minsize;i++){
		data[i]=dft[i][0];
	}
	return minsize;
}

uint DFT::getRealParts(Buffer &b) const{
	uint minsize=min(size,b.size());
	b.clearChannel();
	for(uint i=0;i<minsize;i++){
		b[i]=dft[i][0];
	}
	return minsize;
}

double *DFT::getImaginaryParts() const{
	double *ret=(double*)malloc(sizeof(double)*size);
	for(unsigned int i=0; i<size;i++){
		ret[i]=dft[i][1];
	}
	return ret;
}

uint DFT::getImaginaryParts(double* data, uint dsize) const{
	uint minsize=min(size,dsize);
	for(uint i=0;i<minsize;i++){
		data[i]=dft[i][1];
	}
	return minsize;
}

uint DFT::getImaginaryParts(Buffer &b) const{
	uint minsize=min(size,b.size());
	b.flush();
	for(uint i=0;i<minsize;i++){
		b[i]=dft[i][1];
	}
	return minsize;
}

double *DFT::getMagnetudes() const{
	double *ret=(double*)malloc(sizeof(double)*size);
	for(unsigned int i=0; i<size;i++){
		ret[i]=sqrt(dft[i][0]*dft[i][0]+dft[i][1]*dft[i][1]);
	}
	return ret;
}
uint DFT::getMagnetudes(double* data, uint dsize) const{
	uint minsize=min(size,dsize);
	for(uint i=0;i<minsize;i++){
		data[i]=sqrt(dft[i][0]*dft[i][0]+dft[i][1]*dft[i][1]);
	}
	return minsize;
}

uint DFT::getMagnetudes(Buffer &b) const{
	uint minsize=min(size,b.size());
	b.flush();
	for(uint i=0;i<minsize;i++){
		b[i]=sqrt(dft[i][0]*dft[i][0]+dft[i][1]*dft[i][1]);
	}
	return minsize;
}

void DFT::setData(double* real,double* img, int Size){
	if(dft!=0) free(dft);
	size=Size;
	dft=(fftw_complex*)fftw_malloc(sizeof(fftw_complex)*size);
	for(unsigned int i=0;i<size;i++){
		dft[i][0]=real[i];
		dft[i][1]=img[i];
	}
}

DFT DFT::operator=(const DFT& d2){
	dft=(fftw_complex*)fftw_malloc(sizeof(fftw_complex)*size);
	memcpy(dft,d2.dft, sizeof(fftw_complex)*size);
	return *this;
}
DFT DFT::operator*(const DFT& d2){
	int u = max(size, d2.size);
	DFT ret(u);
	if(size != d2.size){
		
		fftw_complex* res;
		double *re1=getRealParts();
		double *im1=getImaginaryParts();
		double *re2=d2.getRealParts();
		double *im2=d2.getImaginaryParts();
		double *real=resize(re1, size, re2,d2.size);
		double *img=resize(im1, size, im2,d2.size);
		if(size > d2.size){
			res=multiply(re1,im1,real,img,size);
		}else{
			res=multiply(real,img,re2,im2,d2.size);
		}
		fftw_free(ret.dft);
		ret.dft=res;
		free(re1);free(re2);free(im1);free(im2);
		free(real);free(img);
	}else{
		for(unsigned int i=0;i<size;i++){
			ret.dft[i][0]=dft[i][0]*d2.dft[i][0] - dft[i][1]*d2.dft[i][1];
			ret.dft[i][1]=dft[i][1]*d2.dft[i][0] + dft[i][1]*d2.dft[i][0];
		}
	}
	return ret;
}

void DFT::operator*=(const DFT& d2){
	if(size != d2.size){
		double*real,*img;
		fftw_complex* res;
		double *re1=getRealParts();
		double *im1=getImaginaryParts();
		double *re2=d2.getRealParts();
		double *im2=d2.getImaginaryParts();
		real=resize(re1, size, re2,d2.size);
		img=resize(im1, size, im2,d2.size);
		if(size > d2.size){
			res=multiply(re1,im1,real,img,size);
		}else{
			res=multiply(real,img,re2,im2,d2.size);
		}
		fftw_free(dft);
		dft=res;
		free(real);free(img);
		free(re1);free(re2);free(im1);free(im2);
	}else{
		double img=0.0;
		double real=0.0;
		for(unsigned int i=0;i<size;i++){
			real=dft[i][0]*d2.dft[i][0] - dft[i][1]*d2.dft[i][1];
			img=dft[i][1]*d2.dft[i][0] + dft[i][1]*d2.dft[i][0];
			dft[i][0]=real;
			dft[i][1]=img;
		}
	}
}

DFT DFT::operator+(const DFT& d2){
	int u = max(size, d2.size);
	DFT ret(u);
	if(size != d2.size){
		double*real,*img;
		fftw_complex* res;
		double *re1=getRealParts();
		double *im1=getImaginaryParts();
		double *re2=d2.getRealParts();
		double *im2=d2.getImaginaryParts();
		real=resize(re1, size, re2,d2.size);
		img=resize(im1, size, im2,d2.size);
		if(size > d2.size){
			res=add(re1,im1,real,img,size);
		}else{
			res=add(real,img,re2,im2,d2.size);
		}
		ret.dft=res;
		free(re1);free(re2);free(im1);free(im2);
		free(real);free(img);
	}else{
		for(unsigned int i=0;i<size;i++){
			ret.dft[i][0]=dft[i][0]+d2.dft[i][0];
			ret.dft[i][1]=dft[i][1]+d2.dft[i][0];
		}
	}

	return ret;
}

double* DFT::resize(double *part1, unsigned int size1, double* part2, unsigned int size2){
		double delta;
		double *eqx,*real,*eqy;
		int max,min;
		if(size1 > size2){
			max=size1;
			min=size2;
			eqy=part2;
			//interpolate d2
			delta=size1/size2;
			//calculate eqx
			eqx=(double*)malloc(sizeof(double)*size2);
			for(unsigned int i=0;i<size2;i++){
				eqx[i]=i*delta;
			}
			real=(double*)malloc(sizeof(double)*size1);
			
		}else{
			max=size2;
			min=size1;
			eqy=part1;
			//interpolate d
			delta=size2/size1;
			eqx=(double*)malloc(sizeof(double)*size1);
			for(unsigned int i=0;i<size;i++){
				eqx[i]=i*delta;
			}
			real=(double*)malloc(sizeof(double)*size2);
			
		}

		Utility::Interpolate(Utility::MONOTONE_CUBIC, eqx,eqy,min,real, max,1.0);
		free(eqx);
		return real;
}

fftw_complex* DFT::multiply(double* re1,double* im1,double* re2,double* im2,unsigned int size){
	fftw_complex* out=(fftw_complex*)fftw_malloc(sizeof(fftw_complex)*size);
	for(unsigned int i=0;i<size;i++){
		out[i][0]=re1[i]*re2[i] - im1[i]*im2[i];
		out[i][1]=im2[i]*re1[i] + im1[i]*re2[i];
	}
	return out;
}

fftw_complex* DFT::add(double* re1,double* im1,double* re2,double* im2,unsigned int size){
	fftw_complex* out=(fftw_complex*)fftw_malloc(sizeof(fftw_complex)*size);
	for(unsigned int i=0;i<size;i++){
		out[i][0]=re1[i]+re2[i];
		out[i][1]=im2[i]+im1[i];
	}
	return out;
}

void DFT::operator+=(const DFT& d2){
	if(size != d2.size){
		double*real,*img;
		fftw_complex* res;
		double *re1=getRealParts();
		double *im1=getImaginaryParts();
		double *re2=d2.getRealParts();
		double *im2=d2.getImaginaryParts();
		real=resize(re1, size, re2,d2.size);
		img=resize(im1, size, im2,d2.size);
		if(size > d2.size){
			res=add(re1,im1,real,img,size);
		}else{
			res=add(real,img,re2,im2,d2.size);
		}
		fftw_free(dft);
		dft=res;
		free(real);free(img);
		free(re1);free(re2);free(im1);free(im2);
	}else{
		for(unsigned int i=0;i<size;i++){
			dft[i][0]+=d2.dft[i][0];
			dft[i][1]+=d2.dft[i][1];
		}
	}
}

fftw_complex* DFT::getContent() const {
	return dft;
}

DFT::operator Buffer(){
	double *re=getRealParts();
	return Buffer(re,size);
	free(re);
}

}
